library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)
logmod <- function(x) sign(x) * log(1 + abs(x))
sqrtmod <- function(x) sign(x) * sqrt(abs(x))
cbrtmod <- function(x) sign(x) * (abs(x)**(1 / 3))
df <- read.csv("../data/preprocessed_illusion1.csv") |>
mutate(
Block = as.factor(Block),
Illusion_Side = as.factor(Illusion_Side)
)
for (ill in c(
"Ebbinghaus",
"VerticalHorizontal",
"MullerLyer"
)) {
print(ill)
model <- brms::brm(
Error ~ t2(Illusion_Difference, Illusion_Strength, bs = "tp") + (1 | Participant),
data = filter(df, Illusion_Type == ill),
family = "bernoulli",
algorithm = "sampling",
refresh = 400
)
name <- paste0("gam_", tolower(ill), "_err")
assign(name, model) # Rename with string
save(list = name, file = paste0("models/", name, ".Rdata"))
model <- brms::brm(
brms::bf(
RT ~ t2(Illusion_Difference, Illusion_Strength, bs = "tp") + poly(ISI, 2) + (1 | Participant),
sigma ~ t2(Illusion_Difference, Illusion_Strength, bs = "tp") + poly(ISI, 2) + (1 | Participant),
beta ~ t2(Illusion_Difference, Illusion_Strength, bs = "tp") + poly(ISI, 2) + (1 | Participant)
),
data = filter(df, Illusion_Type == ill, Error == 0),
family = "exgaussian",
init = 0,
algorithm = "sampling",
refresh = 400
)
name <- paste0("gam_", tolower(ill), "_rt")
assign(name, model) # Rename with string
save(list = name, file = paste0("models/", name, ".Rdata"))
}
Ebbinghaus
Model Selection
test_models <- function(data) {
# TODO: add random effect
models_err <- list()
models_rt <- list()
for (diff in c(
"Illusion_Difference",
"logmod(Illusion_Difference)",
"sqrtmod(Illusion_Difference)",
"cbrtmod(Illusion_Difference)"
)) {
for (ill in c(
"abs(Illusion_Strength)",
"logmod(abs(Illusion_Strength))",
"sqrtmod(abs(Illusion_Strength))",
"cbrtmod(abs(Illusion_Strength))"
)) {
print(paste(diff, "x", ill))
err <- glmmTMB::glmmTMB(
as.formula(
paste0(
"Error ~ Illusion_Effect / (",
diff,
" * ",
ill,
") + (1|Participant)"
)
),
data = data,
family = "binomial"
)
models_err[[paste(diff, "*", ill)]] <- err
rt <- glmmTMB::glmmTMB(
as.formula(
paste0(
"RT ~ Illusion_Effect / (",
diff,
" * ",
ill,
") + poly(ISI, 2) + (1|Participant)"
)
),
data = data
)
models_rt[[paste(diff, "*", ill)]] <- rt
}
}
mutate(performance::test_performance(models_err), Outcome = "Error") |>
rbind(mutate(performance::test_performance(models_rt), Outcome = "RT")) |>
select(-Model, -log_BF) |>
datawizard::convert_na_to(select = "BF", replacement = 1) |>
arrange(Outcome, desc(BF)) |>
display(footer = "Each model is compared to 'Illusion_Difference'")
}
test_models(filter(df, Illusion_Type == "Ebbinghaus"))
Error Rate
data <- filter(df, Illusion_Type == "Ebbinghaus")
Descriptive
plot_desc_errors_obs <- function(data, n_groups=5) {
data |>
mutate(
Illusion_Difference = datawizard::categorize(Illusion_Difference,
n_groups = n_groups,
split = "equal_length",
label = "median"
),
bins = datawizard::categorize(Illusion_Strength, split = "equal_length", n_groups = 20)
) |>
group_by(Illusion_Difference, bins) |>
summarize(
Illusion_Strength = mean(Illusion_Strength),
Error = sum(Error) / n()
) |>
ggplot(aes(x = Illusion_Strength)) +
geom_bar(aes(y = Error, fill = Illusion_Difference),
alpha = 1 / 3, position = "dodge", stat = "identity", width = diff(range(data$Illusion_Strength)) / 20
)
}
plot_desc_errors <- function(data) {
dat <- data |>
mutate(
Illusion_Difference =
datawizard::categorize(Illusion_Difference,
n_groups = 5,
split = "equal_length",
label = "median"
)
)
col <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(dat$Illusion_Difference)))
plot_desc_errors_obs(data) +
geom_smooth(data=dat,
aes(y = Error, color = Illusion_Difference, fill = Illusion_Difference),
method = "gam",
formula = y ~ s(x, bs = "cr", k=8),
method.args = list(family = "binomial"),
alpha = 1 / 3
) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
scale_color_manual(values = col) +
scale_fill_manual(values = col) +
coord_cartesian(ylim = c(0, 1), xlim = range(dat$Illusion_Strength)) +
labs(x = "Illusion Strength", y = "Probability of Error") +
theme_modern() +
facet_grid(~Illusion_Side, labeller = "label_both")
}
plot_desc_errors(data)
Model Specification
formula <- brms::bf(
Error ~ Illusion_Effect / (logmod(Illusion_Difference) * abs(Illusion_Strength)) +
(1 + Illusion_Effect / (logmod(Illusion_Difference) * abs(Illusion_Strength)) | Participant),
family = "bernoulli",
decomp = "QR"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
illusion1_ebbinghaus_err <- brms::brm(formula,
data = data,
refresh = 100,
normalize = FALSE
)
save(illusion1_ebbinghaus_err, file = "models/illusion1_ebbinghaus_err.Rdata")
Model Inspection
load("models/illusion1_ebbinghaus_err.Rdata")
load("models/gam_ebbinghaus_err.Rdata")
plot_model_errors <- function(data, model, gam) {
pred <- estimate_relation(model, length = c(3, 100))
dat <- data |>
mutate(
Illusion_Difference = datawizard::categorize(
Illusion_Difference,
n_groups = 3,
split = "equal_length",
label = "median"
),
Illusion_Strength = as.character(Illusion_Strength)
) |>
group_by(Illusion_Strength, Illusion_Difference) |>
summarize(
Error = sum(Error),
n = n()
) |>
group_by(Illusion_Strength) |>
mutate(
n = sum(n),
Error = Error / n,
Illusion_Strength = as.numeric(Illusion_Strength)
) |>
ungroup()
plot_desc_errors_obs(data, n_groups = 3) +
scale_fill_manual(values = c("#F44336", "#FFC107", "#4CAF50"), guide = "none") +
ggnewscale::new_scale_fill() +
geom_ribbon(data=pred, aes(
ymin = CI_low, ymax = CI_high,
fill = Illusion_Difference, group = Illusion_Difference
),
alpha = 0.2
) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_hline(yintercept = c(0.5), linetype = "dotted", alpha = 0.5) +
geom_line(
data = estimate_relation(gam, length = c(3, 100)),
aes(y=Predicted, color = Illusion_Difference, group = Illusion_Difference),
linetype = "dotted"
) +
geom_line(data=pred, aes(y=Predicted, color = Illusion_Difference, group = Illusion_Difference)) +
scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
scale_x_continuous(expand = c(0, 0)) +
scale_color_gradientn(colours = c("#4CAF50", "#FFC107", "#F44336"), trans = "reverse") +
scale_fill_gradientn(colours = c("#4CAF50", "#FFC107", "#F44336"), trans = "reverse") +
theme_modern(axis.title.space = 5) +
guides(color = "none") +
labs(
color = "Difficulty",
fill = "Difficulty",
y = "Probability of Error",
x = "Illusion Strength"
) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
p_ebbinghaus_err <- plot_model_errors(data, illusion1_ebbinghaus_err, gam_ebbinghaus_err)
p_ebbinghaus_err
Response Time
data <- filter(df, Illusion_Type == "Ebbinghaus", Error == 0)
Descriptive
plot_desc_rt <- function(data) {
dat <- data |>
mutate(
Illusion_Difference =
datawizard::categorize(Illusion_Difference,
n_groups = 5,
split = "equal_length",
label = "median"
)
)
col <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(dat$Illusion_Difference)))
dat |>
ggplot(aes(x = Illusion_Strength, y = RT)) +
geom_jitter2(aes(color = Illusion_Difference),
width = diff(range(dat$Illusion_Strength)) / 50,
alpha = 1 / 3
) +
geom_smooth(aes(y = RT, color = Illusion_Difference, fill = Illusion_Difference),
method = "gam",
formula = y ~ s(x, bs = "cr", k=8),
alpha = 1 / 10
) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
scale_color_manual(values = col) +
scale_fill_manual(values = col) +
coord_cartesian(ylim = c(0, 2)) +
labs(x = "Illusion Strength", y = "Response Time") +
theme_modern() +
ggside::geom_ysidedensity(aes(fill = Illusion_Difference), color = NA, alpha = 0.3) +
ggside::theme_ggside_void() +
ggside::scale_ysidex_continuous(expand = c(0, 0)) +
ggside::ggside() +
facet_grid(~Illusion_Side, labeller = "label_both")
}
plot_desc_rt(data)
